%%  
% --------------------------------------------------------------
%   1-D, 1-G DIFFUSION-BASED, NONLINEAR RESPONSE MATRIX SOLVER
% --------------------------------------------------------------
%  This test case aims to determine keff for a one-dimensional
%  cartesian slab reactor using analytic response functions that
%  are essentially combinations of transmission/reflection 
%  coefficients and multiples of integrated fluxes.  The current
%  coupling depends implicitly on keff, and a nonlinear problem
%  arises, for which a solution is sought using Newton's method.
%  Additionally, the standard red-black sweeping method has been
%  added for comparison of convergence.  
%  j. roberts, 11/19/2009

clear; clc;
format long

% declare global problem data
global ne dc sa ns sw bcL bcR

problem = 1;

if problem == 1
%------------------------------------------------------------------
% SAMPLE PROBLEM  # 1                                
%  - two slabs of the same material 
%  - zero re-entrant current
% FDM solution yields k ~ 1.002; here, we get 1.0022
ne = 2;                     % number of elements
dc = [ 4.000  4.000 ];      % diffusion coefficient      [cm^2]
sa = [ 0.100  0.100 ];      % absorption cross-section   [1/cm^2]
ns = [ 0.200  0.200 ];      % nu * fission cross-section [1/cm^2]
sw = [ 4.247  4.247 ];      % slab width                 [cm]
% dc = [ 4.000  4.000 4.000  4.000  ];      % diffusion coefficient      [cm^2]
% sa = [ 0.100  0.100 0.100  0.100];      % absorption cross-section   [1/cm^2]
% ns = [ 0.200  0.200 0.200  0.200];      % nu * fission cross-section [1/cm^2]
% sw = [ 2.1235  2.1235 2.1235  2.1235];      % slab width                 [cm]
bcL = 0;                    % left boundary condition
bcR = 0;                    % right boundary condition
k = 1.00;                    % initial guess of keff
j = ones(1,ne*2);           % initial guess for currents
%------------------------------------------------------------------
elseif problem == 2
%------------------------------------------------------------------
% SAMPLE PROBLEM  # 2                                
%  - five element, infinite core, |mod|f1|mod|f2|mod|
%  - refected bc's
% FDM solution yields k ~ 1.0104; here, we get k ~ 1.0102
ne = 5;
dc = [ 0.200  0.650  0.200  0.600  0.200  ];      
sa = [ 0.020  0.140  0.030  0.150  0.020  ];      
ns = [ 0.000  0.175  0.000  0.170  0.000  ];      
sw = [ 2.000  2.500  3.000  2.500  2.000  ];  
bcL = 1;                    
bcR = 1;  
k = 1.;
j = ones(1,ne*2);
%------------------------------------------------------------------
elseif problem == 3
% SAMPLE PROBLEM # 3
%  - twenty element, finite core, something asymmetrical 
ne = 50;
% dc = [ 0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ... 
%        0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ... 
%        0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ... 
%        0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ];
% sa = [ 0.020  0.140  0.220  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020  ...
%        0.020  0.140  0.220  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020  ...
%        0.020  0.140  0.220  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020  ...
%        0.020  0.140  0.220  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020 ];
% ns = [ 0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ...  
%        0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ...
%        0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ...
%        0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ];  
% sw = [ 2.200  2.500  9.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ...  
%        2.200  2.500  9.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ... 
%        2.200  2.500  9.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ... 
%        2.200  2.500  9.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ];
mxel = 140;
tmp = ones(mxel,1); dc = tmp; sa = tmp; ns = tmp; sw = tmp;
% dc(1:2:mxel,1) = 0.24;    dc(2:2:mxel,1) = 0.66;
% sa(1:2:mxel,1) = 0.02;    sa(2:2:mxel,1) = 0.13898687904;
% ns(1:2:mxel,1) = 0.00;    ns(2:2:mxel,1) = 0.15449421059;
% sw(1:2:mxel,1) = 1.25;    sw(2:2:mxel,1) = 1.587;
dc(1:2:mxel,1) = 0.24;    dc(2:2:mxel,1) = 0.66;
sa(1:2:mxel,1) = 0.02;    sa(2:2:mxel,1) = 0.13967685149;
ns(1:2:mxel,1) = 0.00;    ns(2:2:mxel,1) = 0.15449421056;
sw(1:2:mxel,1) = 1.25;    sw(2:2:mxel,1) = 1.587;
% dc(1:2:mxel) = (rand(mxel/2,1)-.5)*.24/2+.24; 
% dc(2:2:mxel) = (rand(mxel/2,1)-.5)*.66/2+.66;
% sa(1:2:mxel) = (rand(mxel/2,1)-.5)*.02/2+.02; 
% sa(2:2:mxel) = (rand(mxel/2,1)-.5)*0.13898687904/2+0.13898687904;
% ns(1:2:mxel) = .0; 
% ns(2:2:mxel) = (rand(mxel/2,1)-.5)*0.15449421059/2+0.15449421059;
% sw(1:2:mxel,1) = 1.25;    sw(2:2:mxel,1) = 1.587;
bcL = 1;                    
bcR = 1;  
k = 1.01;
j = ones(1,ne*2);   
end

%------------------------------------------------------------------
%---------------------- BEGIN SOLVERS -----------------------------
%------------------------------------------------------------------

% normalize initial currents and form vector of unknowns
j       = j/sum(j);
x       = [j k]';
tol     = [1.e-11 1.e-11]; % tolerance for norms


%---------------------- RED BLACK SOLVER --------------------------
% this solver is not set up optimally from a memory standpoint,
% meaning there are lots of zeros that could be eliminated if
% i wanted to be clever.
mxrbit = 20;

tic
[x1,itrb,nrmrb,xx1] = redblack(ne,dc,sa,ns,sw,x,tol,mxrbit);
t1 = toc;
% the relative nonlinear residual, for plotting
plotrb = nrmrb/nrmrb(1);
% the "solution"
xsol = xx1(:,end); 
% estimate rate of convergence
b=1;
p1 = log( norm(xx1(:,end-b-1)-xsol)/ norm(xx1(:,end-b-2)-xsol) ) ...
        / log( norm(xx1(:,end-b-2)-xsol)/ norm(xx1(:,end-b-3)-xsol) );
% compute the final norm
n1=norm( feval('respfct',x1) );    

%---------------------- NONLINEAR SOLVERS -------------------------
% the nonlinear function is   
%   -------------------------------     -----
%  | R(k_n)*M - I  |       0       |   | j_n |   
%  |    F(k_n)     | -L(k_n,j_n-1) | x | k_n | = 0
%   -------------------------------     -----
%  where F is 'gains' and L is 'losses'

%---------------------- (1) NEWTON-GMRES SOLVER -------------------
% uses newton-gmres functions written by kelley

tic
nrm0 = 1;% nrmrb(end);
[x2, ithistGM, ierr, nrmgm, xx2] = nsolgm(x,'respfct',tol,nrm0);
if nrmgm(end) == 0
    nrmgm(end) = eps; % in case the last delta is ~ 0.
end
t2 = toc;
itgm = length(ithistGM);
% the relative nonlinear residual, for plotting
plotgm = ithistGM(:,1)/ithistGM(1,1);
% estimate rate of convergence using RB solution as ref
b = 0;
p2 = log( norm(xx2(:,(end-b-1))-xsol)/ norm(xx2(:,(end-b-2))-xsol) ) ...
        / log( norm(xx2(:,(end-b-2))-xsol)/ norm(xx2(:,(end-b-3))-xsol) );

%---------------------- (2) NEWTON SOLVER -------------------------
% uses the standard newton method w/ a fd-approximated jacobian
x3 = x; xx3(:,1) = x;
tic
z       = respfct(x3); % initial residual
itmx    = 25; 
it      = 1;
zz(1)   = norm(z);
k1      = k; k0 = 0;
nrmcn(1) = 1;
%warning off
%m = mprecond('respfct',x1);
while ( zz(it) > tol(1)*zz(1)+tol(2) && it < itmx )
    % compute F'(x)
    fp          = jacob('respfct',x3,z);
    %disp(['toc 1 = ',num2str(toc)])
    %fp          = diffjac(x3,'respfct',z);
    %disp(['toc 2 = ',num2str(toc)])
    % solve F'*s = -F(x); note fp'fp for spd
    s           = -(fp'*fp)\(fp'*z);
    %disp(['toc 3 = ',num2str(toc)])
    % create new x
    %k0          = x3(end);
    %x0          = x3;
    x3          = x3 + s;
    %k1          = x3(end);
    % normalize currents
    x3(1:end-1) = x3(1:end-1)/sum(x3(1:end-1));
    it          = it+1; 
    xx3(:,it)   = x3; % Keep all x's to estimate rho
    z           = respfct(x3);
    zz(it)      = norm(z); % THIS IS PLOTTED  
end
%warning on
t3 = toc;
% estimate rate of convergence using RB as ref
p3 = log( norm(xx3(:,end-1)-xsol)/ norm(xx3(:,end-2)-xsol) ) ...
        / log( norm(xx3(:,end-2)-xsol)/ norm(xx3(:,end-3)-xsol) );
% relative nonlinear residuals for plotting
plotcn = zz/zz(1);

%--- OUTPUT
semilogy(0:itrb-1,plotrb,'k',0:itgm-1,plotgm,'b--',...
         0:it-1,plotcn,'r-.')
axis([0,max(itgm,max(itrb,it)), ...
    min(min(plotcn(end),plotgm(end)),plotrb(end)), 1.01]);     
xlabel('Iterations')
ylabel('Relative Nonlinear Residuals')
legend('red-black','newt-gmres','newton')

  disp(' ')
  disp(' *** final results ***')
  disp(' ')
  disp('         |    red-black    |    newt-gmres    |    classic      |   ')
fprintf('      it |      %3i        |      %3i         |     %3i         | \n', ...
    itrb-1,itgm-1,it-1)
fprintf('    keff |%16.13f |%16.13f  |%16.13f | \n', x1(end),x2(end),x3(end))

n2=norm( feval('respfct',x2) );
n3=norm( feval('respfct',x3) );
fprintf('||F(x)|| |    %4.3e   |    %4.3e    |    %4.3e   | \n',n1,n2,n3)
fprintf('       p |      %4.3f      |      %4.3f       |      %4.3f      | \n',p1,p2,p3)
fprintf('    time |      %4.3f      |      %4.3f       |      %4.3f      | \n',t1,t2,t3)



% this is the jacobian for the first problem with some analytically
% derived values; the matrix is singular to working precision
% AA = [ -1       0       0.7999247629           0       -0.3917587446e-1*x(2) 
%         0       -1      0.2197291477    0       -0.4031917547e-1*x(2)
%         0       -1      0.2197291477    0       -0.4031917547e-1*x(2)
%         0      0.7999247629  0              -1       -0.3917587446e-1*x(2)
%        -0.1824010047   0.8395989952  0.8395989952  -0.182  -2*(1.894426949*x(1)+0.8944269490*x(2))
%        ]
